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Abstract. Data assimilation methodologies are designed to incorporate noisy observations of 
a physical system into an underlying model in order to infer the properties of the state of the 
system. Filters refer to a class of data assimilation algorithms designed to update the estimation 
of the state in a on-line fashion, as data is acquired sequentially. For linear problems subject to 
Gaussian noise filtering can be performed exactly using the Kalman filter. For nonlinear systems it 
can be approximated in a systematic way by particle filters. However in high dimensions these 
particle filtering methods can break down. Hence, for the large nonlinear systems arising in 
applications such as weather forecasting, various ad hoc filters are used, mostly based on making 
Gaussian approximations. The purpose of this work is to study the properties of these ad hoc 
filters, working in the context of the 2D incompressible Navier-Stokes equation. By working in this 
infinite dimensional setting we provide an analysis which is useful for the understanding of high 
dimensional filtering, and is robust to mesh- refinement. We describe theoretical results showing 
that, in the small observational noise limit, the filters can be tuned to perform accurately in 
tracking the signal itself (filter stability), provided the system is observed in a sufficiently large 
low dimensional space; roughly speaking this space should be large enough to contain the unstable 
modes of the linearized dynamics. Numerical results are given which illustrate the theory. In a 
simplified scenario we also derive, and study numerically, a stochastic PDE which determines filter 
stability in the limit of frequent observations, subject to large observational noise. The positive 
results herein concerning filter stability complement recent numerical studies which demonstrate 
that the ad hoc filters perform poorly in reproducing statistical variation about the true signal. 
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Assimilating large data sets into mathematical models of time-evolving systems presents a major 
challenge in a wide range of applications. Since the data and the model are often uncertain, a 
natural overarching framework for the formulation of such problems is that of Bayesian statistics. 
However, for high dimensional models, investigation of the Bayesian posterior distribution of model 
state given data is not computationally feasible in on-line situations. For this reason various ad 
hoc filters are used. The purpose of this paper is to provide an analysis of such filters. 

The paradigmatic example of data assimilation is weather forecasting: we have many complex 
models to predict the state of the atmosphere, currently involving on the order of O{10^) unknowns, 
but we cannot know the exact state of the system at any one time; we thus have to contend 
with an initial condition which is only known incompletely. This is compensated for by a large 
number, currently on the order of O{10^), partial observations of the atmosphere at subsequent 
times. Filters are widely used to make forecasts which combine the mathematical model of the 
atmosphere and the data to make predictions. Indeed the particular method of data assimilation 
which we study here includes, as a special case, the algorithm commonly known as 3DVAR. This 
method originated in weather forecasting. It was first proposed at the UK Met Office in 1986 [T6] . 
and was developed by the US National Oceanic and Atmospheric Administration (NOAA) soon 
thereafter; see [21]. More details of the implementation of 3DVAR by the UK Met Office can be 
found in [17], and by the European Centre for Medium- Range Weather Forecasts (ECMWF) in 
[6]. The 3DVAR algorithm is prototypical of the many more sophisticated filters which are now 
widely used in practice and it is thus natural to study it. The reader should be aware, however, 
that the development of new filters is a very active area and that the analysis here constitutes an 
initial step towards the analyses required for these more recently developed algorithms. For insight 
into some of these more sophisticated filters see [271 El 1211 [HI ttS] and the references therein. 

Filtering can be performed exactly for linear systems subject to Gaussian noise: the Kalman 
filter [12]. For nonlinear or non-Gaussian scenarios the particle filter [8] may be used and provably 
approximates the desired probability distribution as the number of particles is increased [1]. 
However in practice this method performs poorly in high dimensional systems [23j. Whilst there 
is considerable research activity aimed at overcoming this degeneration [291 [3] , the methodology 
cannot yet be viewed as a provably accurate tool within the context of the high dimensional 
problems arising in geophysical data assimilation. In order to circumvent problems associated 
with the representation of high dimensional probability distributions some form of ad hoc Gaussian 
approximation is typically used to create practical filters, and the 3DVAR method which we analyze 
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here is perhaps the simplest example of this. 

In the paper |15] a wide range of Gaussian approximate filters, including 3DVAR, are evaluated 
by comparing the distributions they produce with a highly accurate (and impractical in realistic 
online scenarios) MCMC simulation of the desired distributions. The conclusion of that work is 
that the Gaussian approximate filters perform well in tracking the mean of the desired distribution, 
but poorly in terms of statistical variability about the mean. In this paper we provide a theoretical 
analysis of the ability of the filters to estimate the mean state accurately. Although filtering is 
widely used in practice, much of the analysis of it, in the context of fluid mechanics, works with 
finite-dimensional dynamical models. Our aim is to work directly with a PDE model relevant in 
fluid mechanics, the Navier-Stokes equation, and thereby confront the high-dimensional nature of 
the problem head-on. Study of the stability of filters for data assimilation has been a developing 
research area over the last few years and the paper [2] contains a finite dimensional theoretical 
result, numerical experiments in a variety of finite and (discretized) infinite dimensional systems 
not covered by the theory, and references to relevant applied literature. Our analysis will build 
in a concrete fashion on the approach in [20] and [13] which were amongst the first to study data 
assimilation directly through PDE analysis, using ideas from the theory of determining modes in 
infinite dimensional dynamical systems. However, in contrast to those papers, we will allow for 
noisy observations in our analysis. Nonetheless the estimates in [13] form an important component 
of our analysis. 

The presentation will be organized as follows: in Section |2] we introduce the Navier-Stokes 
equation as the forward model; in Section |3] we formulate data assimilation as a Bayesian inverse 
problem and show how approximate Gaussian filters can be derived, leading to 3DVAR and 
generalizations; in Section H] we introduce notions of stability and prove the Main Theorems 14.31 
and 14.71 concerning filter stability for sufficiently small observational noise; section H] also contains 
derivation of the stochastic PDE (SPDE) and deterministic PDE which may be used to study 
filter stability for 3DVAR when data is aquired frequently in time and the observational noise 
is large; in Section [5] we present numerical results which corroborate the analysis; and finally, in 
Section |6] we present conclusions. The mathematical tools required to follow the arguments in 
this paper comprise basic ideas from probability on Hilbert space, properties of the Navier-Stokes 
equation, analysis of nonautonomous/random dynamical systems, and properties of stochastic 
PDEs. Together these tools facilitate an analysis which allows mathematical study of filtering 
and, we believe, can be built upon to study other problems arising in the filtering of complex 
systems. 
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In this section we establish a notation for, and describe the properties of, the Navier-Stokes 
equation. This is the forward model which underlies the inverse problem which we study in 
this paper. We consider the 2D Navier-Stokes equation on the torus := [0,L) x [0,L) with 
periodic boundary conditions: 

dtu - uAu + u-Vu + Vp = f for all (x, t) G x (0, oo), 
V-u =0 for all (x,t) e T2 X (0,oo), 

u{x,0) = Uo{x) for all X G T^. 

Here u: x (0, oo) is a time-dependent vector field representing the velocity, 

p: X (0, oo) — 7- M is a time- dependent scalar field representing the pressure, /: — )■ is a 
vector field representing the forcing (which we assume to be time-independent for simplicity), and 
1/ is the viscosity. In numerical simulations (see section E]), we typically represent the solution via 
the vorticity w and stream function (; these are related through u = V"'"C and w = V"*" ■ u, where 
= (^2, -dif. We define 

:= I L— periodic trigonometric polynomials u : [0, L)^ — V ■ u = 0, J u{x) dx = 0^ 

and H as the closure of T-L with respect to the (L^(T^))^ norm. We define P : (L^(T^))^ H to 
be the Leray-Helmholtz orthogonal projector. 

Given k = (^1,^:2)"'", define k~^ := (A;2, — /ci)""". Then an orthonormal basis for H is given by 
i/'feiM^ ^ where 

, , , fc"*" f2'Kik ■ x\ 

for k \ {0}. Thus for -u G we may write 

u= ^ Uk{t)i)k{x) 

k&I?\{0} 

where, since m is a real-valued function, we have the reality constraint U-k = —Uk- We then define 
the projection operators P\ : H ^ H and Qx : H ^ H by 

P\u= ^ Uk{t)iJk{x), Qx = I-P\- 

|27rfcP<AL2 

We let Wx = PxH and = QxH. 

Using the Fourier decomposition of u, we define the fractional Sobolev spaces 

H':=lueH: ^ {A7r^\k\''y\uk\'' < oo\ (1) 
I /fcez2\{o} I 
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with the norm : = (Efc(47^^l'^^>fc^^^^ where s G M. We use the abbreviated notation 1 1 till 
for the norm on H^, and | ■ | for the norm on H = H^. 

Applyng the projection P to the Navier-Stokes equation we may write it as an ODE (ordimary 
differential equation) in H; see [U |25l [22] for details. This ODE takes the form 
dti 

— + uAu + B{u, u) = f, u{0) = uo. (2) 

Here, A = —PA is the Stokes operator, the term B{u, u) = P{u ■ Vu) is the bilinear form found by 
projecting the nonlinear term u ■ Vu into H and finally, with abuse of notation, / is the original 
forcing, projected into H. We note that A is diagonahzed in the basis comprised of the {4'k}k&'2\{o}y 
on H, and the smallest eigenvalue of A is Ai = Att'^/L'^. The following proposition is a classical 
result which implies the existence of a dissipative semigroup for the ODE ([2]). See Theorems 9.5 
and 12.5 in [22j for a concise overview and [25| [26] for further details. 

Proposition 2.1. Assume that uq G and f & H. Then ^ has a unique strong solution on 
t G [0, T] for anyT>0: 

nil 

u G L-((0,T);ifi) n L2((0,T);D(A)), - G L\{Q,T)-H). 

Furthermore the equation has a global attractor A and there is K > such that, if Uq G A, then 
sup,>o||ix(t)f <i^. 

We let {^(^t) : H^}t>o denote the semigroup of solution operators for the equation 

([2]) through t time units. We note that by working with weak solutions, ^'(■, t) can be extended to 
act on larger spaces H^, with s G [0, 1), under the same assumption on /; see Theorem 9.4 in [22] . 
We will, on occasion, use this extension of ^E'(-,t) to act on larger spaces. In particular we use it 
in the following proposition which follows from Lemmas 5.3 and 5.5 in [5]. This proposition is key 
to enabling us to show that the Bayesian inverse problem, which underlies the filtering problem of 
interest, is well-posed. 

Proposition 2.2. Let s G (0, 1] and to > 0. Then, for all t > to, and u,v & H, 

IMuM's <to'C{\f\' + H') 

\\^{U, t) - vl/(^;, < \u-v\, \f\)\u - 

The key properties of the Navier-Stokes equation that drive our analysis of the filters 
are summarized in the following proposition, taken from the paper To this end, define 

^(■) = \1'(-,/;,) for some fixed h > 0. Note that the statement here is closely related to the 
squeezing property [22j of the Navier-Stokes equation, a property employed in a wide range of 
applied contexts. 
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Proposition 2.3. Let u E A and v G H^. There is P = L, z/) > such that 

\\-^{u) - *(tOir < exp{/3h)\\u - vf. (3) 
Now let \\u — v\\ < R and assume that 

X>X* := 

^3 V u J 

where c is a dimensionless positive constant. Then there exists t* = t*{\f\,L,i>,X,R) with the 
property that, for all h G (0,t*], there is •y E (0, 1) such that 

\\Q,{^iu)-^{v))f<y'\\u-vf. (4) 

Proof. The first statement is simply Theorem 3.8 from [T3]. The second statement follows from 
the proof of Theorem 3.9 in the same paper, modified at the end to reflect the fact that, in our 
setting, Px6{0) ^ 0. Note also that the constant A appearing on the right hand side of the lower 
bound for A in the statement of Theorem 3.9 in [13] should be Ai (as the proof that follows in [T3] 
shows) and that use of definition of K (see Theorem 3.6 of that paper) allows rewrite in terms of 
K - indeed the proof in that paper is expressed in terms of K. □ 

3. Inverse Problem: Filtering 

In this section we describe the basic problem of filtering the Navier Stokes equation ([2]): estimating 
properties of the state of the system sequentially from partial, noisy, sequential observations of 
the state. We first set up the inverse problem of interest, in subsection 13.11 Then, in subsection 



13. 2[ we describe the full statistical filtering distribution. We describe the approximate Gaussian 
filters, which are the focus of the remainder of the paper, in subsection 13.31 in subsection 13.41 we 
make a brief remark concerning the extension to completely observed systems. We conclude with 
subsection 13.51 which introduces 3DVAR as an example of the approximate Gaussian filter. 

3.1. Setup 

Throughout the following we write N for the natural numbers {1,2,3,...}, and Z+ := N U {0} for 
the non-negative integers {0,1,2,3,...}. Recall that we have defined \Ef(-) = \I'(-, /i) for some fixed 
h > Q. Our interest is in determining u from noisy observations of P\u. We let X denote for 
any s > 1 and define {uj}j(zz+, Uj G X byQ 

Uj+i := ^(^i). (5) 

I With abuse of notation, subscripts j will indicate times, while subscripts k will denote Fourier coefBcients as 
before in order to avoid confusion. The meaning should also be clear in context. 
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Thus Uj = u{jh) where u solves We now let {^j j^gN be a noise sequence in Wx which perturbs 
the sequence {-Px'^jjjeN to generate the observation sequence {yj}j^n in W\ given by 



exactly. The goal of filtering is to determine information about the state Uj from the data Yj. 
Mathematically, it is natural to formulate this as a Bayesian inverse problem, and this viewpoint 
is described in subsection 13.21 

3.2. Filtering Distribution 

In the statistical formulation of the filtering problem we assume that {uq, Yj) is a random variable on 
(X, Wl), defined by specifying the distributions of Uq (the prior) and = {^i}i<i<j (the observational 
noises) which is assumed to be an i.i.d. sequence. The aim is to find the conditional probability 
distribution on uq given a single realization of Yj, and we denote this conditional distribution by 
F{uo\Yj). The filtering distribution, namely the probability distribution on uj given Yj, denoted 
by Fj{uj\Yj), is then found as the image of P('Uol^) under the map \E'(-; j/i). In finite dimensions 
carrying out this program is a straightforward application of Bayes' theorem. Here we show how 
similar ideas may be applied in the infinite dimensional setting. We make the following assumption. 

Assumption 3.1. The prior distribution on uq is a Gaussian Po('Uo) = N{mo,Co), with the 
property that Fq{H) = 1. The observational noise sequence is an i.i.d sequence in Wx, 

independent of Uq, with C,i distributed according to a Gaussian measure N{0, T) on Wx, with T 
strictly positive on Wx. 

Under this assumption the Bayesian inverse problem has a well-defined solution, and exhibits 
well-posedness with respect to the data, in the Hellinger metric (see for a definition). 

Theorem 3.2. Let Assumption \3.1\ hold. The measure F{uo\Yj) is absolutely continuous with 
respect to Po(^o) with Radon-Nikodym derivative given by 



Vj+i ■■= PxUj+i + ^j+i, j e Z+. (6) 
We let Yj = {yi}i-i, the accumulated data up to time t = jh. We assume that uq is not known 




(7) 



Here 




(8) 



Furthermore, F{uQ\Yj) is Lipschitz in Yj with respect to the Hellinger metric. 
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Proof. Note that W\ is a finite dimensional space. The result can thus be deduced from Corollary 
2.2 and Theorem 2.5 in [5]. To apply Corollary 2.2 it suffices to show that \E'(-; z/i): H is 

continuous for any i G Z since then P\^!{-;ih)\ H — )■ W\ is continuous; the required continuity 
follows from the second item of Proposition 12.21 To apply Theorem 2.5 it suffices to show that 
^{■■,ih):H H is polynomially bounded, since then ^^^^-.H — t- M^, and its Lipschitz constant 
with respect to data Yj, are both polynomially bounded. This polynomial boundedness of ^ follows 
from the first item of Proposition 12.21 □ 

The measure Fj{uj\Yj) is then defined by the push-forward under the semigroup \&(-; j/i) and 
Proposition 12.21 gives the following corollary: 

Corollary 3.3. Let Assumption \3. 1\ hold. Then the sequence of measures {Fj{uj\Yj)}j>o is well- 
defined by 

p,(-|r,) = ^(-;j/^)*n-l^,) 

where P(-|l^) is given in Theorem \3.S[ Furthermore, forv ~ Pj(-|Yj), v G with probability one. 



3.3. Approximate Gaussian Filters; Partial Observations 

The measures P('Uo|y^) and Fj{uj\Yj) determined in Theorem 13.21 and Corollary 13.31 are, in practice, 
very hard to compute by statistical sampling methods. Sequential Monte Carlo methods (SMC) 
can in principle be applied to determine the Fj{uj\Yj), but new ideas are required to extend them 
to problems in high dimensions [23]. The measure F{uo\Yj) has the advantage of being specified 
via density with respect to a Gaussian and it is shown in [15] that, building on this fact, Markov 
chain-Monte Carlo (MCMC) methods can be used successfully to approximate P(no|y,) when the 
dynamics of ([2]) are mildly turbulent; but these methods are very expensive indeed, and do not 
exploit the sequential nature of the problem as j is incremented. Sequential methods are attractive 
for online applications and for this reason various ad hoc approximation methods are used which 
lead to tractable sequential algorithms. A commonly made approximation is to impose a Gaussian 
structure on the filtering distributions so that 

Fj{uj\Yj) ^ N{mj,Cj). (9) 

The key question in designing an approximate Gaussian filter, then, is to find an update rule of 
the form 

(10) 
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Because of the linear form of the observations in together with the fact that the noise is mean 
zero-Gaussian, this update rule is determined directly if we use the Gaussian assumption that the 
distribution of Mj+i given Yj is Gaussian: 

Uj+i\Yj ~ N{mj+i,Cj+i). (11) 

In general, even if the approximation ([9]) is a good one, there is no reason to expect ( ITT]) to be 
a good approximation since the distribution on Uj+i\Yj is the pushforward of Fj{uj\Yj), assumed 
Gaussian, under the nonlinear map \E'(-;/;,) and in general only linear transformations preserve 
Gaussianity. However, in this paper, we will simply impose the approximation ffTTj) with 

m^+i = \l>(mj; h) (12) 

and with Cj+i specified exogenously. This then defines the map (fTOl) . as we now show. 

Assumption 3.4. Assume that Ujj^i\Yj is specified by / fli]) for rrij+i given by (Q^j and covariance 
operator Cj+i on H which is strictly positive on Wx and which commutes with A]^ Assume further 
that yj+i\uj^i is given by ^ where the random variable ^j+i is a mean zero Gaussian in Wx with 
covariance T a strictly positive operator on Wx which commutes with A. 

Note that the Gaussian A^(mj+i, Cj+i) factors as the product of two independent Gaussians on 
Wx and W^, because Cj+i commutes with A; to avoid proliferation of notation, we also denote by 
Cj+i the covariance operator restricted to Wx and W^. The factoring as independent products is 
inherited by the resulting Gaussian approximation for Mj_|_i|Yj_|_i as the following charactertization 
shows. 



Theorem 3.5. Let Assumption 3.4 hold. Then Uj+i\Yj^i is Gaussian on H and factors as the 
product of two indepenent Gaussians on Wx and W^. Denoting the mean and covariance of both 
of these independent Gaussians by fhjj^i and Cj+i respectively we have 



rrij+i 



Qa^(%), on W'^ 

C7+\Pa^(%) + r-iy,+i, on Wx. 



Proof. For economy of notation, let v denote the random variable Uj^i\Yj and y the random 
variable yj+i. Under the stated assumptions, {v,y) is a jointly Gaussian random variable and 
we are interested in finding the conditional distribution of v\y (which is the desired distribution 

§ Note that commuting with A is equivalent to being diagonahzable in the same basis as A, namely the 
{ipk}kez'^\{o}- 
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of Uj+i\Yj^i.) Let h'{dv,dy) denote the Gaussian distribution of {v,y) and let h'o{dv,dy) denote 
the Gaussian distribution specified as the independent product of the Gaussian measures fiQ = 
A^(mj_(.i, Cj+i) and A^(0,r). From the properties of finite dimensional Gaussians it is immediate 
that u has density with respect to uq and that 

dv , ^ ( „ 2 l|-n-i |2 



= exp(^--|r 2(y-p^t;)|2 + _|r \y\ 



(t;) oc exp(^-^|r ^■^{y-Pxv)\ 



duQ 

If n is the desired conditional distribution on v\y, then Lemma 2.3 in [5] gives 

dfi 
dfio 

with constant of proportionality depending only on y. Note that fio factors in the desired fashion 
on Wx and W^, and that the change of measure depends only on Pxv; hence fi factors in the 
same fashion. Furthermore, since the change of meausure depends only on (the finite dimensional) 
Pxv G Wx, completing the square gives the expression for the measure in Wx, and in W^ we have 
/i = /iQ. This completes the proof. □ 

It is demonstrated numerically in [15] that Gaussian approximations of the filtering 
distribution such as the one characterized in the preceding theorem are, in general, not good 
approximations. More precisely, they fail to accurately capture covariance information. However, 
the same numerical experiments reveal that the methodology can perform well in replicating the 
mean, if parameters are chosen correctly, even if it is initially in error. Indeed this accurate tracking 
of the mean is often achieved by means of variance inflation - increasing the model uncertainty, here 
captured in the exogenously imposed Cj, in comparison with the data uncertainty, here captured 
in the F. The purpose of the remainder of the paper is to explain, and illustrate, this phenomenon 
by means of analysis and numerical experiments. To this end we introduce a compact notation for 
the mean update fhj m^+i. 

We define an operator Bj: Wx x W^ ^ Wx x W^ by 

5, = ( ^^^^^^ (13) 

where Cj+i and Cj+i denote the covariances in Wx (resp. W^) in the top left (resp. bottom right) 
entries of Bj. We also extend yj+i from an element of Wx to an element of H in the canonical 
fashion, by defining it to be zero in W^. Then Theorem 13.51 yields the key equation 

fhj+i = Bj^imj; h) + {I- Bj)yj+i, (14) 

which demonstrates that the estimate of the mean at time j ' + 1 is found as an operator-convex 
combination of the true dynamics applied to the estimate of the mean at time j, and the data at 
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time j + 1. Note that, in the case of a hnear dynamical system where h) is a hnear map, the 
matrix / — Bj is the Kalman gain matrix [12j. 

3.4- Approximate Gaussian Filters; Complete Observations 

We will also study the situation where complete observations are made, obtained by taking A — ?■ oo 
in the preceding analyses. The observations are given by 

yj:=Uj + ij, j G Z+ (15) 

where now ?/-,■, G H and the i.i.d. mean zero Gaussian sequence {■CiljeN is defined by a covariance 
operator F on H. 

It is possible to characterize the full filtering distribution in this setting, but rather technical to 
do so; the example of the inverse problem for the heat equation in [24] illustrates the technicalities 
involved. Because of this we concentrate on studying only the equation for the mean update in 
the approximate Gaussian filter. This takes the form f[T^ where, on the whole of if, we have 

B,=C,^,Cjl,, Cjl, = Cjl, + V~\ (16) 
3.5. Example of an Approximate Gaussian Filter: 3DVAR 

The algorithm described in the previous section yields the well-known 3DVAR method, discussed 
in the introduction, when Cj = C for some fixed operator C. To impose commutativity with A, we 
assume that the operators F and C are both fractional powers of the Stokes operator A, in W\ and 
H respectively. Note that fixing Cj = C implies that Cj = C where 

C-i =C'^ + T-\ inWx 

We choose Aq = £ jj] and set C = in H andT = ct'^Aq'^'^ in Wx- Substituting into the 

update formula f|T^ for rhj and defining rj = a/6, a = ( — Bq{vi) = (/ + rfA^°^)~^rfA^ in Wx 
then (|T3l) gives a constant matrix 

Using this we obtain the mean update formula 

m,+i = 5^(m,) + (/ - B)y,+^. (18) 

Notice that for C, F given as above, the algorithm depends only on the three parameters A, a and 
?7, once the constant of proportionality £ in Aq is set. The parameter A measures the size of the 
II The parameter i forms a useful normalizing constant in the numerical experiments of section [5] 
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space in which observations are made; for fixed wavevector k, the parameter r/ is a measure of the 
scale of the uncertainty in observations to uncertainty in the model; and the sign of the parameter 
a determines whether, for fixed rj and asymptotically for large wavevectors, the model is trusted 
more (a > 0) or less (a < 0) than the data. 

In the case A = oo, the case of complete observations where the whole velocity field is noisily 
observed, we again obtain (fT8|) . with B = Bq^t]) = (/ + rj^ A^°')~^iff A^"" in if. The roles of 7] and a 
are the same as in the finite A (partial observations) case. 

The discussion concerning parametric dependence with respect to varying rj shows that, for 
the example of 3DVAR introduced here, and for both A finite and infinite, variance infiation can be 
achieved by decreasing the parameter rj. We will show that variance infiation does indeed improve 
the ability of the filter to track the signal. 

4. Stability 

In this section we develop conditions under which it is possible to prove stability of the 
nonautonomous dynamical system defined by the mean update equation (fl^ . By stability we here 
mean that, when the noise perturbing the observations is 0(e), the mean update will converge to 
an 0{e) neighbourhood of the true signal, even if initially it is an 0{1) distance from the true 
signal. In subsection 14.11 we study the case of partial observations; subsection 14.21 contains the 
(easier) result for the case of complete observations. The third subsection 14.31 shows how our 
results can be applied to the specific instance of the 3DVAR algorithm introduced in subsection 
13.51 for any a G M, provided rj, which is a measure of uncertainty in the data to uncertainty in 
the model, is sufficiently small: this, then, is a result concerning variance inflation. In the final 
subsection 14.41 we also consider 3DVAR in the case of frequent observations and large ?7, deriving 
an SPDE and a PDE, both of which may be used to study filter stability. 

For simplicity, we will assume a "truth" which is on the global attractor, as in [13]. This is 
not necessary, but streamlines the presentation as it gives an automatic uniform in time bound in 
H^. Recall that || • || denotes the norm in if^, and | ■ | the norm in if; similarly we lift || ■ || to 
denote the induced operator norm — > H^. 

4.1. Main Result: Partial Observations 

In this case we will see that it is crucial that the observation space Wx is sufficiently large, i.e. that 
a sufficiently large number of modes are observed. This, combined with the contractivity in the 
high modes encapsulated in Proposition 12.31 from [13], can be used to ensure stability if combined 
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with variance inflation. 

We study filters of the form given in flT^ and make the following assumption on the 
observations {yj}- We note that this assumption is incompatible with the Gaussian assumptions 
used to derive the filter and we return to this point when describing numerical results in section |5l 

Assumption 4.1. Consider a sequence uj = u{jh) where u{t) is a solution of ^ lying on the 
global attractor A. Then, for some X G (Ai, oo), 

Vj = P^^j + 

for some sequence satisfying supj>i < e. 

We make the following assumption about the family {Bj}, and assumed dependence on a 
parameter rj G M^. Recall that the inverse of rj quantifies the amount of variance inflation. 

Assumption 4.2. The family of positive operators {Bj{ri): — ?■ H^}j>i commute with A, satisfy 
supjyi \\Bj{ri)\\ < 1, and supj^i \\I — Bj{ri)\\ < b for some b E M^, uniformly with respect to rj. 
Furthermore, (/ — Bj{ri))Qx = and there is, for all A > Ai, constant c = c(A) > such that 
supj>i \\PxBj{T])\\ < cr/2. 

We now study the asymptotic behaviour of the filter under these assumptions. 



Theorem 4.3. Let Assumptions 4:J_ and\4^hold, choose any mo G B/^i ('u(O), r) and let (A*,t*) 



be as given in Proposition \2.3[ Assume that A > A*. Then for any h G (0,t*] there is t] sufficiently 
small so that the sequence {rhj}j>Q given by [14^ satisfies, for some a G (0, 1), 

l^j — Uj II < a^r + 2be a\ 

=0 



Hence 

limsup ||mj — UjW < e. 

jr— >oo 1 CL 

Proof. Assumption 14.21 shows that y^+i = Px'^{uj)+^j^i. Recall that in (fH]) yj^i has been extended 
to an element of H, by defining it to be zero in W^, and we do the same with ^j+i. Substituting 
the resulting expression for j/j+i in (IT^ we obtain 

m,+i = B,^{m,) + (J - B,)Px^iuj) + (J - 

but since (/ — Bj)Qx = by assumption we have 

m,+i = B.^im,) + (/ - B,)^iu,) + (J - (19) 
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Note also that 

= B,^{uj) + (/ - B,)^{u,). 

Subtracting gives the basic equation for error propagation, namely 

m,+i - = 5,(^(m,) - ^{uj)) + (/ - (20) 

Since A > A* the second item in Proposition 12.31 holds. Fix a G (7, 1) where 7 is defined in 
Proposition 12. 3[ Assume, for the purposes of induction, that 

i-i 

W^j — Uj II < a^r + 2be a\ 

Define R = 2r noting that the inductive hypothesis implies that, for e sufficiently small, 
Wrhj — UjW < r + 26(1 — a)^^e < R. Applying Px to (|2U]) and using gives 

||P,(%+i-n,-+i)|| < ||P,5,||||(M/(%)-M/(«,))ll + lin(/-5.)l|e 

< c{\)r]'^ exp{(3h/2)\\mj - Uj\\ + be. 

Applying Qx to fl2U]) and using (jl]) give^S 

||QA(m,+i-n,+i)|| < ||i?,||||gA(v&(m,)-^(n,))|| + ||gA(/-i?,)lk 

< 'y\\rhj — UjW + be. 

Now note that, for any w G i?"'", ||w|| = (||Paw|P + HQa^IP)^ < UPaW-'H + HQaW-'H- Thus, by adding 
the two previous inequalities, we find that 



1 -Mj+i|| < (c(A)r/2exp(/3V2) +7)ll^j -^11 + 26e. 



'j+ 

Since 7 G (0, 1) and a G (7, 1), we may choose rj sufficiently small so that 

||mj_|_i — Mj+i|| < aWrhj — Uj\\ + 2be. 

and the inductive hypothesis holds with j 1— )■ j + 1. Taking j — > cxd gives the desired result 
concerning the limsup. □ 

Remark 4.4. Note that the proof exploits the fact that Pj\l/(-) induces a contraction within a 
finite ball in H^. This contraction is established by means of the contractivity of Bj in Wx, via 
variance inflation, and the squeezing property ci/^(-) in Wx, for large enough observation space, 
from Proposition \2.3[ 

^ The term be on the right-hand side of the final identity can here be set to zero because (/ — Bj)Q\ = 0; however 
in the analogous proof of Theorem 14.71 it is present and so we retain it for that reason. 
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There are two important conclusions from this theorem. The first is that, even though the 
solution is only observed in the low modes, there is sufficient contraction in the high modes to 
obtain an error in the entire estimated state which is of the same order of magnitude as the error 
in the (low mode only) observations. The second is that this phenomenon occurs even when the 
initial estimate suffers from an 0{1) error. 

4.2. Main Result: Complete Observations 

Here we study filters of tlie form given in ( fT4|) with observations given by ( fT5i) . In this situation the 
whole velocity field is observed and so, intuitively, it should be no harder to obtain stability than in 
the partially observed case. The proof is in fact almost identical to the case of partial observations, 
and so we omit the details. We observe that, although there is no parameter A in the problem 
statement itself, it is introduced in the proof: as in the previous subsection, see Remark I4.4[ the 
key to stability is to obtain contraction in using the squeezing property of the Navier-Stokes 
equation, and contraction in Wx using the properties of the filter to control unstable modes. 
We make the following assumptions: 

Assumption 4.5. Consider a sequence Uj = u{jh) where u{t) is a solution of ^ lying on the 
global attractor A. Then 

Vj = + 

for some sequence satisfying supj^i < e. 

Assumption 4.6. The family of positive operators {Bj{ri): — t- H^}j>i commute with A, satisfy 
\\Bj{ri)\\ < I, and sup.>]^ ||/ — i?j(r]) II < b for some b G , uniformly with respect to rj. 
Furthermore, for all A > Ai, there is a constant c = c(A) > such that supj^i \\P\Bj{ri)\\ < crj^. 

We now study the asymptotic behaviour of the filter under these assumptions. 

Theorem 4.7. Let Assumptions \4-^ and \4-(^ hold and choose any friQ G MHi{u{Q),r). Then for 
any h G (0,t*], with t* given in Proposition \2.3[ there is rj sufficiently small so that the sequence 
{rhj}j>o given by [I4\ ) satisfies, for some a G (0, 1), 



i-i 

< a^r + 2be'^a\ 

i=0 



Hence 

limsup ||mj — UjW < e. 



1 
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Proof. The proof is nearly identical to that of Theorem 14.31 Differences arise only because we 
have not assumed that (/ — Bj)Q\ = 0. This fact arises in two places in Theorem 14.31 The first is 
where we obtain f lT9|) . However in this case we directly obtain f|T9|) since the whole velocity field is 
observed. The second place it arises is already dealt with in the footnote appearing in the proof of 
Theorem 14.31 when estimating the contraction properties in W^; there we indicate that the proof 
is already adjusted to allow for the situation required here. □ 

Remark 4.8. // sup -^j^ < c?7^ then the proof may be simplified considerably as it is not 

necessary to split the space into two parts, Wx and W^. Instead the contraction of Bj can be used 
to control any expansion in "^i-), provided t] is sufficiently small. 

We observe that the key conclusion of the theorem is the stabilization of the algorithm when 
started at distances ofO{l); the asymptotic bound, although of 0(e), has constant which may 
exceed 1 and so the bound may appear worse than that obtained by simply using the observations 
to estimate the signal. In practice, however, we will show that the algorithm gives estimates of the 
state which improve upon the observations. 

4.3. Example of Main Result: 3DVAR 

We demonstrate that the 3DVAR algorithm from subsection 13.51 satisfies Assumptions 14.21 and 14.61 
in the partially and completely observed cases respectively, and hence that the resulting filters will 
locate the true signal, provided rj is sufficiently small. In the next subsection we also consider the 
limit of frequent observations and 77 large, in which case (S)PDEs can be derived to study filter 
stability. 

Satisfaction of Assumptions 14.21 and 14.61 follows from the properties of 

i?o(r/) = (/ + rfAl^rWAl^ I - Boiv) = (/ + v'Al'^yK 
Note that the eigenvalues of Bo{ri) are 

l + 772(4£7r2|fc|2)'"' 

if Aq = iA. Clearly the spectral radius of BqIt]) is less than or equal to one on Wx or H, 
independently of the sign of a. The difference is just that |A;|2 < A/Ai in the former, and \k\ 
is unbounded in the latter. 

First we consider the partially observed situation. We note that Bj = B and is given by (|T7|) : 

i) 



B 
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so that the Kalman gain-hke matrix I — B is given by 

(/ + 7/2^2")-! 



From this it is clear that (/ — B)Q\ = 0. Furthermore, since the spectral radius of -Bo(^) does 
not exceed one, the same is true of B. Hence for the operator norms from into itself we have 

< 1. Similarly, if a < then 6 := = 1, whilst if a > then b = (l + rj^iiXi)^"^ < 1. 

Thus Theorem 14.31 applies. 

In the fully observed case we simply have Bj = B where B = Bo{r]) defined above on H. 
Again < 1 and if a < then \\I - B\\ =6=1, whilst if a > then b = (l + r]'^{i\iY°'^ < 1- 
Thus Theorem 14.71 applies. Note (see Remark 14.81) . that if a < then the proof of that theorem 
could be simplified considerably because ||-B|| < 1 and in fact supj>^ ||-B|| < crj"^. 

Remark 4.9. Recall that Theorem \4 . 7| gives the asymptotic bound Ce on the state estimate where 
C = jz^- When a < and b = 1 we have C > 1; in this case the asymptotic bound exceeds that 
obtained by simply employing the observations. In the case a > 0, we have b < 1 and it is possible 
that C < 1, meaning that the bound may be of direct value to the practitioner. However, regardless 
of the value of C , we emphasize that the value of Theorem \4. 7| is the stabilization from 0{1) initial 
error, and not the value of the constant C in the asymptotic bound. Our numerics will show that, 
in practice, the error in the state estimator often falls well below the average error committed by 
simply using the data as an estimator. 

4.4. (S)PDE limit 

Theorems I4.3l and l4.7l are valid for sufficiently small rj and hence exploit variance inflation. In this 
limit the observations are given large weight at low wavevectors and the contraction resulting from 
this acts to stabilize the system. There is an interesting family of (S)PDEs for the mean which can 
be derived in the case of large rj., by considering the limit of frequent observations. We describe 
these (S)PDEs in the case of the 3DVAR method from subsection 13.51 For simplicity we consider 
the fully observed case; partial observations can be handled similarly. 

To this end we assume that /i ^ 1, r G (0, 1], and that = cTq/Zi'" and 5^ = uoaQh^^^ . Thus 
Tj^'^ = uh. We then obtain 

m,+i = (/ + r]^Al")-'r]^Al»^{m,) + (/ + 7]^Al")-'y,^i 

= (/ + 7^-^A,^^)-'^{mj) + (/ + 7^-2^0 2")-i77-2Ao-2°y,+i 
= (/ + w/iAo + (/ + w/iAo ^°)-^w/iAo 
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If we define the sequence {zj}j^^+ by 

Zj+i = Zj + hyj+i, Zo = (21) 

then the preceding expression for fhj+i can be rearranged and expanded formally in powers of h 
to give 

Note also that formal expansion in h gives 

h) = u- h{uAu + B{u, u)-f) + 0{h^). 

Thus substituting in the previous expression and taking the limit /i — )■ we obtain the continuous 
time filter 

— + uAm + B{m,m) +uA^'^''(m- —] = f, m(0) = mo- (22) 

Here the observation y = ^ enters as a forcing to the underlying Navier-Stokes equation. The 
term proportional to m — ^ acts to drive the solution towards the observation, and may compete 
with destabilization induced by the Navier-Stokes forcing uArh + B{m,m) — f. As such it is a 
clear continuous time analogue of the filter update equation (IT^ . 

As in the discrete case (HM . we can study the stability of the continuous time filter by 
expressing the noise in terms of an underlying signal u which the filter aims to uncover, and 
study the difference m — u. We assume that u itself solves the Navier-Stokes equation ([2]). We 
now express the observation signal z in terms of the truth u in order to facilitate study of filter 
stability. If, instead of making the assumption of bounded noise as in Assumption 14. 5[ we assume 
a Gaussian noise consistent with derivation of the filter, then we obtain 



where {Awj}j>i is an i.i.d. sequence and Awi ~ A^(0, /) has the distribution of a white noise in 
H. If r = 1 then we have an Euler-Maruyama discretization of the SDE (stochastic differential 
equation) 

— =u + aoA,^^, ^(0) = 0, (23) 

where is a Brownian motion in time with covariance the identity in H. If r G (0, 1) the limit is 
simply the ODE 

^=u, z{0) = 0. (24) 
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If r = 1 then substituting the expression for z from fl25]) into fl22p we obtain the following 
SDE: 

— h z^Am + i3(m, m) + ^"(m — m) = / + wao^Q - — , m(0) = mg- (25) 

On the other hand, if r G (0,1), then substituting the expression for z from ( IMl) into (1221) we 
obtain the following ODE: 

— h z^Am + i3(m, m) + wAq ^"(m — n) = /, m(0) = mo, (26) 

Equations (1251) and fl2B]) can be used to study filter stability in this frequent observations limit. 
The equations are equivalent to an SPDE and PDE of Navier-Stokes type, driven by noise in the 
former case, and with an additional damping term driving the solution towards the signal u in 
both cases. We expect that for large enough u, which is a form of variance inflation in this high 
frequency data context, the solution will indeed stabilize towards the signal u, although the nature 
of white noise forcing means that excursions from the signal in the former case will occur infinitely 
often. These excursions will then be stabilized by the signal. It would be interesting to analyze the 
properties of the SPDE and PDE by using the theory of nonautonomous and random dynamical 
systems, and the ergodicity theory developed in [TO] . 



5. Numerical Results 

In this section we describe a number of numerical results designed to illustrate the range of filter 



stability phenomena studied in the previous sections. We start, in subsection 15.11 by describing 
two useful bounds on the error committed by filters; we will use these guides in the subsequent 
numerics. Subsection 15.21 describes the common setup for all the subsequent numerical results 
shown. Subsection 15.31 describes these results in the case of complete observations in discrete time, 
whilst Subsection 15.41 extends to the case of partial observations, also in discrete time. Subsection 
15.51 studies filter stability in the case of continuous time observations, using the (S)PDEs derived 
at the end of section |H 

Our theoretical results have been derived under Assumptions 14. ll and 14.51 on the errors. These 
are incompatible with the assumption, underpinning derivation of the approximate filters, that 
the observational noise sequence is Gaussian. This is because i.i.d Gaussian sequences will not 
have finite supremum. In order to test the robustness of our theory we will conduct numerical 
experiments with Gaussian noise sequences. 
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5.1. Useful Error Bounds 

We describe two useful bounds on the error which help to guide and evaluate the numerical 
simulations. To derive these bounds we assume that the observational noise sequence is i.i.d 
with E^j = and E^j ®^j = T. Then 

E\^^\' = tT{T) = Y,9k 

k 

where {gk} are the eigenvalues of the operator T. (This operator must be trace class if the Gaussian 
measure A^(0, F), used to derive the approximate filters, is to defined on H). 

• The lower bound is derived from (!20|) . Using the assumed independence of the sequence we 
see that 

E|%+i - > E|(/ - 5,)0+ir = tr((/ - B,)T{I - 5,)*) (27) 

• The upper bound on the filter error is found by noting that a trivial filter is obtained by 
simply using the observation sequence as the filter mean; this corresponds to setting Bj = 
in fll4p . For this filter we obtain 

E\mj+i - nj+ip < E|^j+ip = tr(F) (28) 
in the case of complete observations, and 

E|m,+i - < E|e,+ip + |Qa%+i|' = tr(F) + IQa^.+iI' (29) 

in the case of incomplete observations. 

Although the lower bound (!27|) does not hold pathwise, only on average, it provides a useful 
guide for our pathwise experiments. The upper bounds (1251) and fl2^ do not apply to any numerical 
experiment conducted with non-zero Bj, but also serve as a useful guide to those experiments: it 
is clearly undesirable to greatly exceed the error committed by simply trusting the data. We will 
hence plot the lower and upper bounds as useful comparitors for the actual error incurred in our 
numerical experiments below. We note that, for the 3DVAR example from subsection 13.51 with 
complete observations, the upper and lower bounds coincide in the limit — )■ as then B ^ 0. 
For partial observations they differ by the second term in the upper bound. 

5.2. Experimental Setup 

For all the results shown we choose a box side of length L = 2. The forcing in Eq. ([2]) is taken 
to be / = V'^ip, where ip = cos{7Tk ■ x) and V"*" = JV with J the canonical skew-symmetric 
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matrix, and k = (5,5). The method used to approximate the forward model ([2D is a modification 
of a fourth-order Runge-Kutta method, ETD4RK |7] , in which the Stokes semi-group is computed 
exactly by working in the incompressible Fourier basis {'ipk{x)}kez'^\{o}j ^"^^ Duhamel's principle 
(variation of constants formula) is used to incorporate the nonlinear term. Spatially, a Galerkin 
spectral method [11] is used, in the same basis, and the convolutions arising from products in the 
nonlinear term are computed via FFTs. We use a double-sized domain in each dimension, buffered 
with zeros, resulting in 64^ grid-point FFTs, and only half the modes in each direction are retained 
when transforming back into spectral space in order to prevent aliasing, which is avoided as long 
as fewer than 2/3 of the modes are retained. 

The dimension of the attractor is determined by the viscosity parameter u. For the particular 
forcing used there is an explicit steady state for all z/ > and for u > 0.035 this solution is stable 
(see [18], Chapter 2 for details). As u decrease the flow becomes increasing complex and the regime 
u < 0.016 corresponds to strongly chaotic dynamics with attributes of turbulent scalings in the 
spectrum. We focus subsequent studies of the filter on a mildly turbulent (z/ = 0.01) parametric 
regime. For this small viscosity parameter, we use a time-step of 6t = 0.005. 

The data is generated by computing a true signal solving (|2]) at the desired value of z/, and 
then adding Gaussian random noise to it at each observation time. Such noise does not satisfy 
Assumption 14.51 since the supremum of the norm of the noise sequence is not finite, and so this 
setting provides a severe test beyond what is predicted by the theory; nonetheless, it should be 
noted that Gaussian random variables only obtain arbitrarily large values arbitrarily rarely. 

All experiments are conducted using the 3DVAR setup and it is useful to reread the end of 
subsection 13.51 in order to interpret the parameters a and rj. We consider both the choices a = ±1 
for 3DVAR, noting that in the case a = —1 the operator B has norm strictly less than one and so 
we expect the algorithm to be more robust in this case (see Remark 14.81 for discussion of this fact). 
For all experiments we set i = X^^ which ensures that the action of Aq", and hence B, on the 
first eigenfunction is independent of the value of a; this is a useful normalization when comparing 
computations with a = 1 and a = — 1. 

In the discrete time experiments we set the observational noise to white noise F = a^I (i.e. 
/9 = in section 13. 5p . Here a = 0.04, which gives a standard deviation of approximately 10% 
of the maximum standard deviation of the turbulent dynamics. Since we are computing in a 
truncated finite-dimensional basis the eigenvalues are summable; the situation can be considered 
as an approximation of an operator whose eigenvalues decay rapidly outside the basis in which we 
compute. 
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5.3. Complete Observations; Discrete Time 

We start by considering discrete and complete observations and illustrate Theorem 14 .7^ and in 
particular the role of the parameter rj. The experiments presented employ a large observation 
increment oi h = 0.5 = 1005t. For a = 1 we find that when rj = a (Fig. [1]) the estimator stabilizes 
from an initial 0{1) error and then remains stable. The upper and lower bounds are satisfied (the 
upper bound after an initial rapid transient), and even the high modes, which are slaved to the low 
modes, synchronize to the true signal. For t] = 10a (Fig. |2]) the estimator fails to satisfy the upper 
bound, but remains stable over a long time horizon; there is now significant error in the k = (7, 7) 
mode, in contrast to the situation with smaller rj shown in Fig. [T] Finally, when rj = lOOcx (Fig. 
|3]), the estimator really diverges from the signal, although still remains bounded. 

When a = —1 the lower and upper bounds are almost indistinguishable and, for all values of 
r] examined, the error either exceeds or fiuctuates around the upper bound; see Figures SJ E] and 
|6] where rj = a, 10a and lOOo" respectively. It is not until r] = 100a (Fig. [6]) that the estimator 
really loses the signal. Notice also that the high modes of the estimator always follow the noisy 
observations and this could be undesirable. For both rj = 100a and lOa, the a = — 1 estimator 
performs better than the one for a = 1 in terms of overall error, illustrating the robustness alluded 
to in Remark 14.81 since for a < we have ||i?|| < 1. However, an appropriately tuned a = 1 filter 
has the potential to perform remarkably well, both in terms of overall error and individual error 
of all modes (see Fig. [1], in contrast to Fig. H]). In particular, this filter has an expected error 
substantially smaller than the upper bound, which does not happen for the case of a = — 1 when 
complete observations are assimilated. 

5.4. Partial Observations; Discrete Time 

We now proceed to examine the case of partial observations, illustrating Theorem 14. 3[ Note that 
the forced mode has magnitude = 50, so ensuring that it is observed requires that A > 50Ai. 
When enough modes are retained, for example when A = lOOAi in our setting, the results for the 
a = 1 case remain roughly the same and are not shown. However, in the case a = — 1, in which the 
observations are trusted more than the model at high wavevectors, the results are greatly improved 
by ignoring the observations of the high-frequencies. See Fig. [71 This improvement, and indeed 
the improvement beyond setting B = for both cases a = ±1 disappears as A is decreased. In 
particular, when A = 25Ai the error is never very much smaller than the upper bound. This is due 
to the fact that the dynamics of the low wavevectors tend to be unpredictable and they contain 
very little useful information for the assimilation. Then, for much smaller A = 4Ai, once enough 
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unstable modes are left unobserved, there is no convergence. The order of magnitude of the error 
in the asymptotic regime as a function of rj remains roughly consistent as A is decreased until the 
estimator no longer converges. For small h (high-frequency in time observations) and complete 
observations, the estimator can be slow to converge to the asymptotic regime. In this case, the 
number of iterations until convergence, for a sufficiently small rj, becomes significantly larger as A 
is decreased (again until the estimator fails to converge at all). 

Given A ~ ^a-^i, we expect that for rj sufficiently small the contribution of the model to the 
filter will be negligible for all k with \k\ < kx for a = 1. Hence the estimators for both a = ±1 will 
behave similarly. An example of this is shown in Fig. |8]where rj = O.Ola = 0.0004 and A = 49Ai in 
Fig. [HI In both cases, the estimator is essentially utilizing all the available observations. There are 
enough observations to draw the higher wavevectors of the estimator closer to the truth than if we 
just set the population of those modes to zero. In contrast, as mentioned above, when A = 25Ai, 
there are not enough observations even when they are all used, and the error is roughly the same 
as the upper bound as r/ — )• (not shown). 

5.5. Continuous Observations 

Finally, we explore the case of continuous observations using the SPDE and PDE derived in 
subsection 14.41 We let a = 1/2 throughout. We invoke a split-step scheme to solve Eqns. [25] and 
[26] in which we compose solutions of ([2]) and the Ornstein-Uhlenbeck process 

^ + wAg ^"(m -u) = tjo-o^o "^(0) = ^0, (30) 

at each step, when r = 1, and setting ctq = in (!30]l when r < 1. We begin by considering the 
r = 1 case in which we recover the SPDE (p5l) . Notice that the parameter u sets a time-scale for 
relaxation towards the true signal, and o"o sets a scale for the size of fiucuations about the true 
signal. The parameter /3 rescales the fiuctuation size at different wavevectors with respect to the 
relaxation time. First we consider setting /3 = 0. In Fig. [9] we show numerical experiments with 
u = 100 and ctq = 0.05. We see that the noise level on top of the signal in the low modes is almost 
0(1), and that the high modes do not synchronize at all; the total error remains 0(1) although 
trends in the signal are followed. On the other hand, for the smaller value of ctq = 0.005, still with 
u = 100, the noise level on the signal in the low modes is moderate, the high modes synchronize 
sufficiently well, and the total error is small. See Fig. [TO] 

Now we consider the case /3 = 1. Again we take u = 100 and ctq = 0.05 and 0.005 in Figures 
[TT] and [121 respectively. The synchronization is stronger than that observed for /3 = in each case. 
This is because the noise decays more rapidly for large wavevectors when /3 is increased, as can be 
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observed in the relatively smooth trajectories of the high modes of the estimator. 

For the case when r < 1 and we recover a PDE, the values of ctq and /3 are irrelevant. The value 
of u is the critical parameter in this case. For values of u of O(IOO) the convergence is exponentially 
fast to machine precision. For values of a; of 0(1) the estimator does not exhibit stable behviour. 
For intermediate values, the estimator may approach the signal and remain bounded and still an 
0(1) distance away(see the case a; = 10 in Fig. [T3|l . or else it may come close to synchronizing 
(see the case a; = 30 in Fig. [T^ . 

6. Conclusion 

This paper contains three main components: 

• we show that the filtering problem for the Navier-Stokes equation may be formulated as a 
well-posed inverse problem in function space, and demonstrated how various approximate 
Gaussian filters can be derived (Theorems 13.21 and 13. 5[ and Corollary 13. 3p : 

• we prove Theorems 14.31 and 14. 7^ which establish filter stability for small enough observational 
noise; 

• we derive an SPDE (l25l) and a PDE ( l26l) . which may be used as the basis to study filter 
stability in the case of frequent observations and large observational noise. 

Numerical results are used to illustrate, and extend the validity of, the theory. 
There are a number of natural future directions which stem from this work: 

• to develop analogous filter stability theorems for more sophisticated filters, such as the 
extended and ensemble-based methods, when applied to the Navier-Stokes equation; 

• to rigorously derive, and study the properties of, the SPDE (1251) which characterizes filter 
stability in the case of frequent observations and large observational noise; 

• to study model-data mismatch by looking at filter stability for data generated by forward 
models which differ from those used to construct the filter; 

• to study the effect of filtering in the presence of model error, by similar methods, to understand 
how this may be used to overcome problems arising from model-data mismatch. 
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Figure 1. Example of a stable trajectory for 3DVAR with v — 0.01, h ~ 0.5, -q ^ a = 0.04, a = 1. 
The top left plot shows the norm-squared error between the estimated mean, m{tn) — m„, and 
the signal, u{tn), in comparison to the preferred upper bound (i.e. the total observation error 
tr(r) = S) and the lower bound tr[(/ — Bn)T{I — The other three plots show the estimator, 

m{t), together with the signal, u{t), and the observations, yn for a few individual modes. 
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Figure 2. Example of a destabilized trajectory for 3DVAR with the same parameters as in Fig 
except the larger value of 77 = lOo" = 0.4. Panels are the same. 
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Figure 3. Example of a destabilized trajectory for 3DVAR with the same parameters as in Fig 
except the larger value of ?7 = lOOcr = 4. Panels are the same. 
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Figure 4. Example of a stable trajectory for 3DVAR with the same parameters as in Fig. [l]except 
with a = — 1. Panels are the same. 
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Figure 5. Example of a stable trajectory for 3DVAR with the same parameters as in Fig. 
with value of a = — 1. Panels are the same. 



FIGURES 



32 




Figure 6. Example of a destabilized trajectory for 3DVAR with the same parameters as in Fig. |3] 
except with value of a = — 1. Panels are the same. 
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Figure 7. Example of an improved estimator for partial observations with A = lOOAi and otherwise 
the same parameters as in Fig. [Sj Panels are the same. 
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Figure 8. Examples of estimators for partial observations with A = 49 and -q — Q.la — 0.004, 
otherwise the same parameters as in Figs. [T]and|31 Left panels are for a = 1 and right panels are 
for a = — 1. 




Figure 9. Trajectories of various modes of the estimator fh and the signal u are depicted above 
for (3^0 and (Tq = 0.05, along with the total relative error in the P norm, \fh — 
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Figure 10. Trajectories of various modes of the estimator m and the signal u are depicted above 
for /3 = and do = 0.005, along with the relative error. 
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Figure 11. Trajectories of various modes of the estiraator ffi and the signal u are depicted above 
for /3 = 1 and cto = 0.05, along with the relative error. 
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Figure 12. Trajectories of various modes of the estimator ffi and the signal u are depicted above 
for /3 = 1 and cto = 0.005, along with the relative error. 
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Figure 13. Trajectories of various modes of the estimator m and the signal u are depicted above 
for r < 1 and cj = 10, along with the relative error. 




Figure 14. Trajectories of various modes of the estimator fh and the signal u are depicted above 
for r < 1 and uj = 30, along with the relative error. 



